##########################################################################################

library(data.table)
library(parallel)
library(Seurat)
library(Rmagic)
library(optparse)
library(reticulate)
use_python("~/miniconda3/envs/magic_v1/bin/python", required = NULL)

##########################################################################################

if(1!=1){
    ## 若未按照magic
    library(reticulate)
    use_python("~/miniconda3/envs/magic_v1/bin/python", required = NULL)
    #wget https://cran.r-project.org/src/contrib/Archive/Rmagic/Rmagic_2.0.3.tar.gz
    install.packages("~/tools/Rmagic_2.0.3.tar.gz")
}

##########################################################################################
option_list <- list(
    make_option(c("--input_dir"), type = "character"),
    make_option(c("--out_path"), type = "character") 
)

if(1!=1){

    ## 输出
    out_path <- "~/20231121_singleMuti/results/seurat_raw"

    input_dir <- "~/20231121_singleMuti/results/cellranger_arc"

}

###########################################################################################
parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

input_dir <- opt$input_dir
out_path <- opt$out_path

dir.create(out_path , recursive = T)

###########################################################################################
## 文件路径
file1_path <- paste0(input_dir , "/testis01/outs/filtered_feature_bc_matrix")
file2_path <- paste0(input_dir , "/testis02/outs/filtered_feature_bc_matrix")
file3_path <- paste0(input_dir , "/testis03/outs/filtered_feature_bc_matrix")

###########################################################################################
## 创建seruat对象
### testis01 ###
all_matrix_01 <- Read10X(data.dir = file1_path)
exp_matrix_01 <- all_matrix_01[["Gene Expression"]]  ## 36601 gene, 20000 cell ##
rna_01 <- CreateSeuratObject(counts = exp_matrix_01, assay = 'RNA', min.cells = 0, min.features = 0, project = '10x_RNA')  

### testis02 ###
all_matrix_02 <- Read10X(data.dir = file2_path)
exp_matrix_02 <- all_matrix_02[["Gene Expression"]]  ## 36601 gene, 9815 cell ##
rna_02 <- CreateSeuratObject(counts = exp_matrix_02, assay = 'RNA', min.cells = 0, min.features = 0, project = '10x_RNA')  

### testis03 ###
all_matrix_03 <- Read10X(data.dir = file3_path)
exp_matrix_03 <- all_matrix_03[["Gene Expression"]]  ## 36601 gene, 11725 cell ##
rna_03 <- CreateSeuratObject(counts = exp_matrix_03, assay = 'RNA', min.cells = 0, min.features = 0, project = '10x_RNA')  

###########################################################################################
### 数据合并 ###
rna_all <- merge(rna_01, y = c(rna_02, rna_03), add.cell.ids = c("testis01", "testis02", "testis03"), project = "testis")  ## 41540 cell ##
mito.features <- grep(pattern = "^MT-", x = rownames(x = rna_all), value = TRUE)
## 计算每个细胞中线粒体基因的表达量占总表达量的百分比。
percent.mito <- Matrix::colSums(x = GetAssayData(object = JoinLayers(rna_all) , slot = 'counts')[mito.features, ]) / Matrix::colSums(x = GetAssayData(object = JoinLayers(rna_all), slot = 'counts'))
rna_all$percent.mito <- percent.mito
rna_all@meta.data$cell <- rownames(rna_all@meta.data)
rna_all@meta.data$batch <- sub("(\\w+)_(\\w+)-(\\w+)","\\1",rna_all@meta.data$cell)

###########################################################################################
## 原始的关于质控参数的图
out_file <- paste0( out_path , "/quality_rawcell.pdf" )
pdf(file=out_file,width=8,height=6)
VlnPlot(object = rna_all, features = c("nFeature_RNA", "nCount_RNA", "percent.mito"), group.by="batch",pt.size = 0)
dev.off()

###########################################################################################
## 细胞质控
rna_all_QC3 <- subset(x = rna_all, subset = nFeature_RNA > 500 & nFeature_RNA < 10000 & nCount_RNA > 500 & nCount_RNA < 50000 & percent.mito < 0.25)  

## 对经过质量控制的RNA测序数据 rna_all_QC3 进行归一化处理，以消除不同样本之间的技术差异
rna_all_QC3_magic <- NormalizeData(object = rna_all_QC3)
## 对归一化后的数据进行尺度缩放，使得每个基因的表达值在样本间具有相似的方差
rna_all_QC3_magic <- ScaleData(object = rna_all_QC3_magic)
## 应用 MAGIC（Model-based Analysis of Single-cell Genomics）算法，对数据进行批次效应的去除和缺失值的填补
rna_all_QC3_magic <- magic(JoinLayers(rna_all_QC3_magic), knn=10, t=6, npca=20)

## 根据批次信息将处理后的数据集拆分为多个子数据集，每个子数据集对应一个批次。
testis.list <- SplitObject(rna_all_QC3_magic, split.by = "batch")

## 对每个子数据集进行独立的归一化和变异特征的识别
testis.list <- lapply(X = testis.list, FUN = function(x) {
    x <- NormalizeData(x)
    x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 3000)
})

## 选择在不同批次数据集中反复变化的特征，以进行数据整合
testis_features <- SelectIntegrationFeatures(object.list = testis.list)

## 找到用于整合数据的锚点特征
testis_anchors <- FindIntegrationAnchors(object.list = testis.list, anchor.features = testis_features)

## 利用锚点特征对数据进行整合，生成整合后的数据集
testis_combined <- IntegrateData(anchorset = testis_anchors)

## 将整合后的数据集的默认测定方法设置为 "integrated"
DefaultAssay(testis_combined) <- "integrated"

###########################################################################################
#### 数据预处理
## 对数据进行尺度缩放
testis_combined <- ScaleData(testis_combined, verbose = FALSE)
## 运行主成分分析（PCA）
testis_combined <- RunPCA(testis_combined, npcs = 100, verbose = FALSE)
## 在PCA空间中运行UMAP降维算法
testis_combined <- RunUMAP(testis_combined, reduction = "pca", dims = 1:30)

## 在PCA空间中找到相邻细胞
testis_combined <- FindNeighbors(testis_combined, reduction = "pca", dims = 1:30)
## 基于找到的相邻细胞，根据指定的分辨率进行聚类
testis_combined <- FindClusters(testis_combined, resolution = 0.3)
## 在聚类后的空间中再次运行UMAP算法，以便于可视化聚类结果
testis_combined <- RunUMAP(object=testis_combined, dims = 1:30)

## 画聚类图 ##
DefaultAssay(testis_combined) <- "MAGIC_RNA"

out_file <- paste0( out_path , "/QC3_magic3_cluster.pdf" )
pdf(file=out_file,width=8,heigh=6)
DimPlot(testis_combined, reduction = "umap", label = TRUE)
dev.off()

## marker基因 ##
out_file <- paste0( out_path , "/QC3_magic3_cluster_marker.pdf" )
pdf(file=out_file,width=8,heigh=6)
FeaturePlot(testis_combined, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP","CD8A"))  ## 这个基因名不是睾丸的marker基因,你要画的话需要改一下
dev.off()

###########################################################################################
## 保存seruat对象
out_file <- paste0( out_path , '/testis_combined.Rdata' )
save( testis_combined, file = out_file )
